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ABSTRACT 


Continuous-time models of oscillator phase noise xCt) usually have 
stationary nth differences, for some n. The covariance structure of such a 
model can be characterized In the time domain by the structure function: 
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* E x(s+t) a” x(s) 
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Although formulas for the special case D2(0 ;t, t) (the Allan variance times 
2 t^) exist for power-law spectral models, certain estimation problems require 
a more complete knowledge of (0) . 

We exhibit a much simpler function of one time variable, D(t), from 
which (0) can easily be obtained from the spectral density by uncomplicated 
Integrations. Believing that D(t) Is the simplest function of time that 
holds the same information as (0), we call D(t) the fundamental structure 
function. 

We compute D(t) for several power- law spectral models. Two examples are 
D(t) = Klt|3 for random walk FM, D(t) = Kt^ Inltj for flicker FM. Then, to 
demonstrate Its use, we exhibit a BASIC program that computes means and vari- 
ances of two Allan variance estimators, one of which Incorporates a method of 
frequency drift estimation and removal. Except for a one-line function defi- 
nition of D(t), the program Is Independent of the phase noise spectrum. The 
outputs were used for assigning confidence Intervals to the results of recent 
hydrogen maser performance tests at JPL. 

THE STRUCTURE FUNCTION 

The purpose of this paper Is to demonstrate an easy way of computing a 
class of time-domain oscillator stability measures known collectively as the 
structure function . Let the phase of an oscillator with nominal frequency vq 
be modeled by 


2TrvQ(t + xCt)) , 

where x(t) is a random process representing the "phase time" of the 


*The research described in this paper was carried out by the Jet Propulsion 
Laboratory, California Institute of Technology, under contract with the 
National Aeronautics and Space Administration. 
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oscillator. The most restricted version of the structure function, discussed 
by Lindsey and Chie [1], is 

D^(T) = E[A^ f , (1) 

where is the backwards difference operator defined by A^f (t) = f (t) - f (t-t) , 
and E denotes the mathematical expectation. It is to be understood that the 
nth difference process A“ x(t) is wide-sense stationary, so that (1) does not 
depend on s. For n = 2 we obtain the Allan variance as 

^ D,(x) . 


A more general version of the structure function is 

Tj^,T 2) = E A^ x(s+t) a“ x(s) , (2) 

1 2 

which was used by Yaglom [2] as a basis for a theory of processes with 
stationary nth differences. For the application to be given at the end of 
this paper, we shall need the yet more general version 

D(t; a,b) * E A^x(s+t) A|^x(s) , (3) 

where 

JL ~ » • * • » a^) , A^ = Ag^ ... A^ , 

— 1 n 

and likewise for A^^. 

By using the transfer function 1 - e of the operator Ax, one can 
express all of these quantities as iiitegrals Involving the two-sided spectral 
density of the process x(t) , which is asstuaed to have stationary nth 

differences for some n.. For example, if the long-term frequency drift rate is 
zero, then 


-/”|1 - S ^ O .) 

— oo 

00 

= 16 / sin^(|-WT) s^(to) ll , 

~oo 

D(t; a,b,c,d) = E A^A^^ x(s+t) A^A^ x(s) 

= fa - e-^“®3(l - e-^“b(l - - e^“^) e^“t S^(m) || . 


(4) 


( 5 ) 


(This is (3) with n » 2.) 
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Because of the Importance of the Allan variance as a stability measure, 
the integral (4) has been evaluated in closed form for power-law spectral 
models 


SjjCw) = a * -2, -1, 0, 1, 2, 

with a high-frequency cutoff when a ^ 1 ("PM” noises) [3], The expressions 
(5) and (3) are more complex, yet we shall show how to compute any of these in 
two easy steps; 1) By means of simple integrations involving only the spec- 
tral density, evaluate a single function iD(t) of one titoe variable; 2) apply 
a certain difference operator of order 2n to D(t) . This function D(t), 
depending only on the spectral density, will be called the fundamental struc- 
ture function because, merely by taking linear combinations of values of D(t), 
one can generate the entire covariance structure of the differences of x(t). 

STATIONARY PHASE 

If x(t) is stationary, then D(t) turns out to be just the autocovariance 
function; 

oo 

D(t) = Cov(x(s+t),x(s)) » f e^“*^ 

*.00 

Indeed, by expanding the differences in the middle expression of (5), one 
can show directly that 

D(t; a,b,C,d) - \ D(t) . (7) 

It is just this form of expression that will be extended to the nonstationary 
situation, where S (u) is not integrable near zero frequency. 

THE FUNDAMENTAL STRUCTURE FUNCTION IN GENERAL 

Let us be given a process x(t) with stationary nth differences. It is 
known [2] that x has a spectral density S (u) with the properties 

X 

/ S„(w) do) < «» , 

1 


f 0^ S„(oj) dw < » , (8) 

0 * 

for some Integer k, 0-k ^ 2n. (We ignore the more general case of a spectral 
distribution function.) Here is one way to compute the fundamental structure 
function D(t) ; Define the function 

B(z) = / e^“* (im)^S^(m) || (9) 
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on the upper half-plane Im z > 0. Let C(z) be any kth integral of B(z), i.e., 
any function such that (z) = B(z) on Im z > 0. It will always be found 
that C(z) can be extended to the real line by continuity. Then let 

D(t) = 2 Re C(t), t real (10) 

If k > 0, then D(t) is not unique; any pol 3 momial of degree < k may be added 
to it. If X is stationary then k = 0, and D(t) reduces to the covariance 
function (6). 

RANDOM WALK FM 

To see how easy it is to carry out the above procedure for power-law 
noise types, consider the spectral density and take k = 4. Then 

BU, 

a fourth integral of which is 

(One can throw out a polynomial of degree less than 4.) Taking In z * 
ln[z| + i Arg z for Im z Z 0 gives 

D(t) = 2 Re C(t) = 0 (t > 0) 

- -t^/6 (t < 0). 

3 

Adding t /12 to this gives the alternate form 

D(t) = It 1^/12 . 


OTHER POWER-LAW NOISE TYPES 

We now give D(t) for the power-law spectrum 


S (oj) = K 1 01 i ^ , K “ - — — — 
^ “ « 2(2ir)“ 


The forms for the constant and the power of oi link the results to an 
accepted notation for the noise spectrum, namely 

■ V“ 

[3], where y = dx/dt, f = oi/(2it), and Sy(f) is the one-sided spectral density 
of y. Fractional values of a are allowed. We can take k to be the integer 
part of 2 - a (refer to (8)). 
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Case 1. 


a < 1, not an odd integer. 


D(t) = K 


- t 


il-a 


a 2r(2— a) cos(Tra/2) 

This case includes white FM (a = 0) and random walk FM (a = -2). 
Case 2. a = -1, -3, -5,... 

K. 


( 11 ) 


^ ' TT ^ (1-a) ! 


(. 12 ) 


This case, which includes flicker FM (a = -1) , is actually the easiest of all. 
Case 3. a = 1 (Flicker PM) . 

Here we use the spectrum 

1-1 


S (m) = exp (- loj/m^l) . 


The exponential high-frequency cutoff yields the elementary result 

K 

D(t) = - -^ ln(t^ + l/o)^) , (13) 

whereas the usual rectangular cutoff yields a cosine integral. 

Case 4. a > 1 . 

This phase noise process (provided with a high-frequency cutoff) is 
stationary; thus D(t) is just the autocovariance function of x. 

DERIVING THE STRUCTURE FUNCTION FROM D(t) 

Having the fundamental structure function D(t) in hand, we now show how 
to use it. Let x(t) be a noise with spectral density Sx(to), k the smallest 
integer satisfying (8) , and ng the smallest integer such that 2ng k k. The 
differences of x of order ng and higher are stationary, and to simplify 
matters let us assume that they are ergodic and have mean zero. For ng = 2, 
this means that an oscillator with phase time x(t) has no long-term average 
frequency drift. Our main result is the following formula for the general 
structure function (3) of order n Z ng; 


D(t; a,b) = A^A_^D(t). (14) 

The proof, which has appeared elsewhere [4], starts from the spectral 
integral 
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(15) 


J ^OO 

^a^-b 
00 — — 
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in which the differences operate on el“t as a function oft. (Safer to (5) 
for the case n » 2.) One cannot simply pull the difference operators outside 
the integral, for the Fourier transform integral of diverges (unless 

k « 0) if the singularity at ai « 0 is not neutralized. One way to do this is 
through the function B(z) given in (9). Another embodiment of D(t) is given 
directly by the formula 


D(t) 



00 

[cos a>t 


1 


1 + ( 0 ‘ 


,2n 


E 


j-0 


(2j)I ^ 




( 1 «) 


where n ^ n^ . 

MOMENTS OF STABILITY ESTIMATORS 

The remainder of this paper presents a nontrivial application of the 
preceding theory. Let the phase time x(t) have a long^tni^ frequency drift 
component : 

x(t) ^ XQ(t) + -Ict^ , 

where c is the constant rate of frequency drift, and XQ(t) is a Gaussian 
process whose 2nd differences are stationary, ergodic, and have mean zero. 
Observing x(t) for 0 ^ t ^ T, we wish to estimate c and the quantities 

2 12 2 

0 (t) “ — E [A X (t)] (Gross Allan variance) , 

^ 2x 

o^.(t) » E [A^ x«(t)]^ (Net Allan variance) , 
yu ^ ^ 

which are the theoretical Allan variances before and after removal of drift. 

To define and manipulate the estimators, it is necessary to set up some 
notation as follows: Write 

C(a,b,t) “ A^A^x(t) . 

(This has nothing to do with (10).) Then EC ® c for any a,b,t. A particular 
one of these is used for estimating c, namely 

c - C(t^,T-t^,T) (17) 
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where Tq = T/6.29, a value chosen to minimize the variance of c in the 
presence of flicker FM noise [5]. The interpretation of (17) is that the 
estimated drift rate equals the average frequency near the end of the record, 
minus the average frequency near the beginning, divided by the intervening 
t ime T - X j, . 

Let us introduce the additional quantities 
c. = C(T,T,jx), c^ = 

where j is an integer. To estimate (2 /t 2) 
statistic 


m 

j=2 


C(t,T-t,T) , 


(18) 


a2(x), we use the unbiased 




( 19 ) 


where it is now assumed that T/x is an Integer m. This is just the usual 
Allan variance estimator involving the sum of squares of adjacent second 
differences. (The unusual scale factor 2/x2 is convenient for these calcu- 
lations.) To estimate (2/x^) a^Q(x), we use 


m 


j=2 


(c, - ty 


( 20 ) 


Since 


we have 


m 


j=2 


V - 2cc^ + c 

T 


( 21 ) 


Our goal is to compute the mean and variance of Vq, a biased estimator. 
First, we see that Vq does not depend on the true drift rate c, and so for 
this purpose we can take c - 0. Then 


EV = i2/xh <J^q(x) 


( 22 ) 
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Var V 


1 




(23) 


(m-1) ‘ 


L 2^ C°v(cj .=k) 


j k 


EVq = EV - 2Ecc^ + Ec , 


(24) 


Var Vq - Var V + 4 Var(cc^) + Var - 4 Cov(V,Sc^) + 2 Cov(V,S^) 


-4 Cov(fic^,S^) . 


(25) 


Because x(t) is Gaussian, the fourth moments in (25) are determined from 
the second moments of c, Cj, c^, and the general formula 

Cov(uv,xy) ■ Eux Evy + Euy Evx , (26) 

where u, v, x, y are zero-mean random variables with a joint Gaussian 
distribution. In turn, the required second moments are just special cases 
of the structure function formula (14), which here takes the form 

E C(a,b,s + t) C(c,d,s) ■ 

where D(t) is the fundamental structure function of the phase‘-tlme noise 
x(t) . This expression, when written out, is a sum of 16 terms involving D. 
(See line 1040 of the BASIC program in Fig. 1.) 

A BASIC PROGRAM 

The tedious but elementary task of computing (23) (25) by means of 

(26) and (27) is carried out by the BASIC program shown in Fig. 1. The main 
point to observe about this program is that the whole algorithm is Independent 
of the phase noise spectrum, which enters only through the one-line function 
definition of D(t) in line 1030 (and the print statement 130). The only 
requirement on S 3 j(io) is that the integer k in (8) be at most 4. This guaran- 
tees that x(t) has stationary second differences. Because of the way the out- 
puts are scaled, a multiplicative constant in D(t) can be neglected; thus we 
can let D(t) = |t|^ for random walk FM. For flicker FM, line 1030 becomes 

DEF DD(T) = T*T*LOG(ABS(T) + (T=0)) . 

Line 1040 evaluates (27), and the subsequent definitions evaluate the required 
special cases. Note that the length of the test run is used as a unit of 
time. It must be admitted that computational efficiency has been sacrified to 
gain ease of coding. 

Figure 2 shows the output of the program for random walk FM. The column 
MEAN (NET) gives the expected value of (t2/2)Vq relative to the true net Allan 
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variance c'yoC'^) • The variances of V and Vg, the estimators of gross and net 
Allan variances, are presented in terms of "degrees of freedom" (DF), defined 
for a positive random variable X by 


DF-iieai . 

Var X ’ 

as if X had a distribution. The DF(GIWSS) result is valid only if the true 
drift has been subtracted from the phase data. When T/t » 1, the gross and 
net DF are both 1, since both V and Vq are squares of mean-zero Gausslans. 

CONFIDENCE INTERVALS FOR ALLAN VARIANCES 

Following Howe et al. [6], we use the results shown in Fig. 2 to assign 
confidence intervals based on the distribution, usually with a fractional 
DF. The author has not estimated the error caused by pretending that V and Vq 
are proportional to x^ variables. Random walk FM was used as a noise model for 
T > 10^ s. Fig. 3 shows the estimates of gross- Allan variance (before drift 
removal) and net Allan variance (after drift removal) for a 72-day test run of 
a pair of hydrogen masers at JPL. 

Different drift estimators are used in Figs. 3a and 3b. In Fig. 3a, we 
use a value of drift measured by retuning the masers over a period much longer 
than the 72-day test run. Regarding this value as close to the "true" value 
c, and using it in (21) in place of c, we can use the gross DF numbers from 
Fig. 2 to assign^confldence intervals to net Allan variance. In Fig. 3b, the 
estimated value c was used, so that one must use the net mean and DF numbers 
to compute the confidence Intervals. The negative bias of Vq pushes the con- 
fidence intervals up; in fact, for T/t 2 (one sample of second phase differ- 
ence) , the bias is so great that the 90% confidence Interval does not contain 
the estimate Itself. This is pushing things too far, perhaps. 

The basic problem is that one cannot remove estimated drift (via (21), 
for example) without also taking a bite out of the long-term random fluctua- 
tions. The outputs of the above computer program for phase noise 

(-2.5 S a ^ 0) support the conjecture that the bias of the net Allan variance 
estimator Vq is always negative. 
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iOO ! MOMENTS OF ALLAN MARXANCE ESTIMATORS 
110 GOSUB 1000 
120 MCW6.29 

130 PRINT N'NOISE TYPE » RANDOM WALK FM' 

140 PRINT 'MEAN t* D.F. FOR GROSS «, NET A.M. '\\ 

' T/T AU S ' ME AN < NET ) ' , ' D . F . ( GROSS ) ' D . F . < NET > ' 

ISO FOR M«2 TO 9,10 TO 18 BY 2,20 TO SO BY 5 
160 GOSUB 1080 

170 DF«2/VARV,DFO«23ICEVOHCEVO/VARVO 
180 PRINT N,EVO,DF,DF0 
190 NEXT M 
200 STOP 

1000 ! DEFINITIONS 

1030 DEF DD<T)»ABS<T#T#T) IFUNDAMENTAL STRUCTURE FUNCTION 

!FOR RANDOM WALK FM 

1040 DEF MM < A , B , C , D , T >» < DD < T ) ~DD < T~A > ~DD < T~B )-DD < T+C > -DD < T+D > 

+DD < T~A~B > +DD < T~ A+C ) +DD < T-A+D > +DD < T-B+C ) +DD< T-B+D ) +DD <T+C+D ) 
-DD < T-A-B+C ) -DD ( T-A-B+D ) -DD < T-A+C+D > -DD < T-B+C+D ) 

•♦•DD < T-A-B-fC-i-D > ) / < A#BJKC*D > 

1050 DEF FF(A,B)»MM<A,1-A,B,1-B,0) 

1060 DEF GG<A,B,J>»MM(A,l-A,B,B,i-j5llB) 

1070 DEF HH(A,J)»MM<A,A,A,A,J3KA> 

1073 RETURN 

1080 ! COMPUTATION SUBROUTINE 

! INPUTS; M«T/TAU, MC«T/TAUC <«6,29) 
lOUTPUTS: EVO/A.V., VAR<V>/A.V.#*2, VAR<V0)/A.V.»:*2 
EVO VARU VAR VO 

1090 U«i/M,Ut>i/MC !U STANDS FOR TAU 
1100 FUU»FF<U,U> ,FCU«FF<UC,U> ,FCC«FF<UC,UC) 

1110 HUO»<HH(U,0>,HU02»HUO«HUO !HU0 IS SCALED ALLAN VARIANCE 
1120 ! 

1130 EVO«HUO-2«FCU<^-FCC 

1140 Ml»M-i,X-0\ H*HH<U,J>,X«X-KMi-J>#H#H FOR J«i TO M~2\ 

X»2#X-fM UKHUO 2 , V AR V*XJ|{2/ < M 1 »M i ) 

1150 VARCU«FCC)|cFUU-fFCU#FCU, VARCC*2!ltFCC3|tFCC 

1160 X«0\ X*X-^GG<UC,U,J>#GG<U,U,J> FOR J»2 TO M\ COVVCU«2*X/Mi 
1170 X«0\ B«GG<UC,U,J>,X«X-i-G3|{G FOR J«2 TO M\ C0VVCC«2*X/Mi 
1180 CVCUCCaRJicFCCJitFCU 
1190 ! 

1200 VARV0aVARV+4*VARCU^fVARCC-4J|tCOMVCU-^2J|tCOVVCC-4#CVCUCC 
1210 RETURN WHERE EV0»EV0/HU0 , VARV«VARM/HU02, VARVO«VARV0/HUO2 


Figure 1. A BASIC program for computing mean and variance 
of an estimator of net Allan variance. 
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NOISE TYPE « RANDOM WALK FM 
MEAN & D.F. FOR GROSS 4. NET A.V. 


T/TAU 

MEAN (NET) 

D.F. (GROSS) 

D.F. (NET) 

2 

11213718 

1 

1. 000 0011 

3 

.4131003 

1 . 882353 

1.2011257 

4 

.56608639 

2 . 7692308 

1.9797428 

S 

.65837896 

3.6571431 

2. 8213698 

6 

. 72007427 

4 . 5454549 

3 . 6927653 

7 

.76417726 

5 . 4339623 

4.5779951 

8 

.7970189 

6.3225806 

5.4662905 

9 

.82222714 

7.2112679 

6.3534235 

10 

.84209356 

8.1000005 

7 . 2390502 

12 

.87125838 

9 .8775517 

9.0083684 

14 

. 89153524 

11.655173 

10.777728 

16 

.90639572 

13.432836 

12.546251 

18 

. 91772997 

15.210527 

14.314574 

20 

.92664775 

16.988236 

16.084209 

25 

. 9423454 

21 . 432559 

20.511747 

30 

.95254386 

25.876923 

24.943548 

35 

.9596919 

30.321313 

29.378236 

40 

.96497606 

34.765708 

33.814985 

45 

.96903914 

39.210128 

38.253179 

50 

.97225997 

43 . 654528 

42.692561 


Figure 2. Output of the program of Fig. 1 for 
random walk FM noise. 
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QUESTIONS AND ANSWERS 


None for Paper #12 
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